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00 \ We propose a strategy for conducting lattice QCD simulations at fixed volume but variable quark 

mass so as to investigate the physical effects of dynamical fermions. We present details of techniques 
which enable this to be carried out effectively, namely the tuning in bare parameter space and efficient 
stochastic estimation of the fermion determinant. Preliminary results and tests of the method are 
presented. We discuss further possible applications of these techniques. 
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Q\ • I. INTRODUCTION 



First results from full simulations of lattice QCD have confirmed the magnitude of the computational task ahead 
and have shown only glimpses of physics beyond the quenched approximation. A recent survey of results can be found 

Preliminary results by the UKQCD collaboration using an 0(a) improved action have shown a surprisingly 
r* ■ strong dependence of the effective lattice volume on the bare quark mass. This complicates chiral extrapolations of 
. ' simulation measurements and obscures comparisons with quenched calculations. With this in mind, we investigate 
. how one might control the effective lattice volume by tuning the bare action parameters while the effects of decreasing 
!-h ' quark mass are studied. 

Before proceeding, we should clarify what we mean by 'effective lattice volume'. For dehniteness, consider the 
Wilson discretisation of QCD giving a lattice action dependent on two bare parameters (3 and K defined in the usual 
way. In the quenched approximation, we have become used to thinking of (3 as uniquely controlling the lattice spacing 
a. At fixed /?, one makes lattice measurements of the rho mass or the Sommer scale ro, defined by |Q 

rA =1.65, 

which is conceptually simpler for the present discussion since, in this case, there is no need to consider extrapolation 
in the (valence) quark mass. One then matches the lattice value of ro to its physical value (« 0.49 fm as extracted 
from heavy quark spectroscopy ||) to obtain the lattice spacing at that value of j3. This mapping olq((3) between j3 
and physical lattice spacing is model dependent in that it is unique to the quenched approximation. The continuum 
limit is not accessible directly since the lattice volume vanishes. One must remain at lattice spacing small enough 
that discretisation errors are small but not so small that finite volume effects are significant. 

There are two ways of extending these ideas in the presence of dynamical fermions controlled by the bare mass 
parameter k. 

1. The conventional procedure is to construct a similar mapping between (3 and lattice spacing a where the 
matching is made using the lattice value of ro extrapolated in k (sea quark mass) to the chiral limit. This yields 
a unique, but regularisation-dependent, mapping a x ([3). Comparisons with the continuum limit are made as in 
the quenched case. 
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2. Alternatively, one may consider matching the lattice value of ro at finite values of k. Here, the picture is that 
the simulation is being done with sea quarks of non-infinite mass. Each value of the bare quark mass (or k) then 
corresponds to a different approximation to continuum QCD with light dynamical quarks, in much the same 
way as does quenched QCD (infinite k). Matching in this case results in a mapping a(f3, n). Clearly, this is also 
regularisation-dependent . 

The term 'effective lattice volume' refers to this second definition of lattice spacing. Our proposal, then, is to conduct 
simulations in which one attempts to hold the effective lattice spacing fixed. In this way one is better able to keep 
the physics constant, control lattice artifacts and finite volume effects while studying the effects of light dynamical 
quarks. In contrast, when adopting the first strategy, the significance of lattice artifacts and finite volume effects 
changes as the chiral extrapolation is made. 

In order to carry out this programme one requires a practical way of identifying curves in the (3, k plane of constant 
lattice spacing 

a(f3, k) — constant. (I) 

Consider F, the lattice measurement of observable / with dimension d, so that 

F = fa~ d . (2) 

Then, in the scaling region and to leading order, curves of constant F yield estimates of these curves of constant 
lattice spacing. In this way, one can track changes in (3 required to compensate for changes in k. As a simple example, 
one can identify the (3 shift involved in comparisons of quenched and dynamical simulations. In practice, there will be 
residual dependence on the choice of F . We would expect ro to be a 'good' choice for exposing sea quark dependence 
whereas m n would not, due to the strong dependence on valence quark mass and the effects of chiral symmetry 
constraints. 

In the rest of the paper, we show how the operator and action matching technology introduced in Q can be used 
to identify such curves. We demonstrate efficient algorithms for achieving it and present some numerical tests. In 
section [n] we summarise the relevant matching formalism required and how it may be used. In section ^ we describe 
an efficient algorithm for making stochastic estimates of the fermion determinant . Results of numerical tests are 
presented in the next section. This is followed by a discussion of additional applications of these techniques, including 
parallel tempering simulations with dynamical fermions. Conclusions and outlook are contained in the final section. 



II. CURVES OF CONSTANT PHYSICS 

We first review the matching formalism introduced in ||. Consider actions Si[U] and 5*2 [U] describing two lattice 
gauge theories with the same gauge configuration space {£/} so that (i = 1, 2) 

Z, = J VUe~ s > [u] , < F >i= ^- J VUFe~ s \ (3) 

For example, Si might be the quenched Wilson action and S2 the 0(a)-improved action for 2-flavour QCD |^|. In the 
present application, we will consider Si and S2 to be the same improved fermion action but at different points in the 
[3, k plane. Here, F is some lattice observable. Expectation values with respect to the two actions can be related via 
a cumulant expansion whose leading behaviour implies j^] 

< F > 2 =< F >i + < FA12 >i +. .. (4) 
where A 12 = Si - S 2 , F = F- < F > etc. (5) 

In general, an action is a function of several parameters. For example, the Wilson action depends on the bare 
parameters {3 and n. In j^] we considered matching action parameters in one of three distinct ways 

Ml: match a given set of operators, i.e. require 

< F n >i=< F n >2 ; 
M2: minimise the 'distance' between the actions, i.e. cr 2 (Ai2); 
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M3: maximise the acceptance in an exact algorithm for 52 constructed via accept/reject applied to configurations 
generated with action Si. 

It was shown that, to lowest order, tuning prescriptions M2 and M3 coincide. In fact, if the operators F n contribute 
to the action with weights which are considered as tuning parameters, then prescription Ml also coincides to lowest 
order. The prescriptions differ in a calculable way at next order. Details are in 
In the present application we take 

Si = S e s(/3 ,K ) and S 2 = 5eff(A k) (6) 

and seek to explore the bare parameter dependence of the lattice theory using configurations generated at a series of 
reference points (/3q, Ko) m parameter space. 

Here, S e f{ is the effective action corresponding to lattice QCD with the fermions integrated out. i. e. 

S cS = -0Wn - T (7) 

where Wu is the usual Wilson plaquette action 

Wn = \ Yl ReTr ^ (8) 
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and 



T = n f TrlnM[U} = ^-Trln^M) . (9) 



The fermion matrix M for the non-perturbative 0{a) improved theory is a function of both k and (3. This is because 
the parameter c sw ^ is a function of (3 |fTo|| . Since the improvement scheme fixes c sw (/3), one must not treat c sw as 
an independently tunable parameter. However, as we shall see, the fact that the operator T(/3, k) is a function of [3 
as well as k introduces some practical complications. 
According to eqn. ^, one requires measurements of 

A12 = S 1 -S 2 = ((3- po)W n + T(J3, k) - T{(3 , k ) (10) 



in order to carry out parameter tuning. We discuss efficient algorithms for this in section III 



For now, consider matching lattice observable F at two neighbouring points in the (3, k plane: 

(A),«o) and {P,k) = ((3 +5(3,k q + Sk). (11) 

According to prescription Ml above and (Q) we require, to first order in small quantities, 

<FA 12 >!=0. (12) 

From this we can deduce that the constant F curve is given by 

S_P = < F6f >r 
5k 5k < FWa >i 

where 

Sf = f(J3 + 8p,K + 5K)-f(Po,K ). (14) 

Equation ( |l3] ) amounts to a non-linear differential equation since the right hand side involves 8(3 (via the c sw param- 
eter). Linearising and taking the limit yields for the constant F curve, 



(15) 



dK d<F> < p { w n + ^ 6sw ) > 



The quantity 



dc sw 
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is well determined [[tO| and so the determination of constant F curves reduces to measuring correlations of the form 

< FW a >i and < FST >i . (16) 



The details of this will be described in section [V 



As pointed out in the Introduction, the details of these curves will depend on the choice of F. For sensible choices 
and reasonably physical values of the parameters, one would hope that the corresponding curves of constant a, o,f(P, k) 



say, would agree rather closely, locally at least. For demonstration purposes, we will consider in section IV several 
simple choices for F: 

• P, the average plaquette, proportional to Wa- 

• Wl, various Wilson loops. 

• Seff, the complete effective action itself. 

• Correlation matrices for measuring the static potential and r$. 

• Hadron correlators. 

The first of these, P, may be readily measured with high accuracy and so is excellent for testing the basic technology. 
However, it is not expected to shed much light on lattice spacing. The last two are computationally more demanding 
but more relevant to the project at hand, identifying curves of fixed physical volume. 

One might expect that matching the full action (F = S e s) would be more physically relevant than matching the 
plaquette piece of the action. From ( |l5| ) we see that this curve is determined by correlations of the form 

< TSf >! (17) 

in addition to those of (Jl6|). As we shall see in section IV, it is more difficult to obtain unbiased estimators for these. 

Now consider matching scheme M2 where 'distances' in the action space are minimised. One can think of this as 
defining 'geodesies' in the (3, k plane with respect to the metric implied by (0) i.e. 

g MV =< (d^idj) > . 

The corresponding affine connection would be 

r MQ/3 =< (d^id^s) > . 

Simple minimisation yields, to first order, 

dd _ <(^ + &<w)f> 



dK < (W n + ^c s 



(18) 



In section ^ we show that these curves are directly relevant to simulations of full QCD using parallel tempering. 
Again, we note that these curves involve operator correlations which are more complex to estimate. 

Finally in this section, we observe that some of the above formalism simplifies considerably in the case of the 
unimproved Wilson action (c sw = 0) since then, for example, 

df 

— =< FWu > . (19) 



III. STOCHASTIC ESTIMATOR OF THE FERMION DETERMINANT 

We require an unbiased estimator for T = TrLni? where H = M^M is a hermitian positive-definite matrix. We 
will also require estimators for T 2 , ST and T5T. Bai, Fahey and Golub 0] have recently proposed estimators, with 
bounds, for quantities of the form 

u f g(H)v (20) 
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where g is some matrix function. In our application g is the logarithm and, for convenience of notation, we set 



L = Ln(H) . (21) 
Taking u — v — 0j, some normalised noise vector (e.g. Zi or Gaussian), we can obtain a stochastic estimate of T via 

E T = —Y,^Mi. (22) 

The corresponding variance of this estimator is 

a\Er) = ^Tr(L 2 ) (23) 

for complex Gaussian noise, and something less than this for Zi noise (±1 on each of the complex component). In 
the case of Gaussian noise, we also obtain rather directly an efficient unbiased estimator for T 2 , 

E T 2 = (E T f - -±-E Q (24) 
where Eq is an unbiased estimator for Q = TrL 2 , 

^ = F^^' (25) 

™ i— 1 

In the case of Zi noise, the corresponding estimator is not readily accessible via the techniques described below, so 
we restrict the discussion to complex Gaussian noise. 

In a companion paper we give fuller details of methods for evaluating the quantities 4>\L(j)i, (f>\L 2 (j)i and so 
on. In practice, we make no use of the bounds presented in §|. Instead we use large enough Lanczos systems so 
that the numerical convergence renders the bounds irrelevant. The efficiency of the method results from an elegant 
relationship between the nodes and weights required for a iV-point Gaussian quadrature and the eigenvalues and 
eigenvectors (so-called Ritz pairs) of a Lanczos matrix of dimension N H . In the companion paper JO]] , we show that 
this relationship, and resulting accuracy remain good even when orthogonality is lost. This is an important point for 
our application. If the Lanczos system is large enough to avoid truncation errors, one is well into the regime where 
orthogonality is lost in standard numerical Lanczos methods. In the present paper we merely summarise the formulae 
required to obtain the present set of results. Preliminary results using these techniques were presented in W . 

The actual estimator which we use for T = TrL is 

AT* 

Et = W Y. 1 ^) (26) 

where 

N 

J(^) = E w I Ln (A J -). (27) 

3=1 

Here {A* } (j = 1,2... N) are the eigenvalues of the TV-dimcnsional tridiagonal Lanczos matrix formed using <pi as 
a starting vector. The weights {i^ 2 } are related to the corresponding eigenvectors In fact, u>j is just the first 
component of the jth eigenvalue of the tridiagonal matrix. 

The estimator Eq, for Q = TrL 2 , is obtained from (§|) and (||) using ln(A 2 ) rather than bi(Xj). For T 2 (see (|J)), 
we define 

E T 2 = (E T ) 2 - -t-tEq (28) 

It is straightforward to show that, with the above definitions 

(£ T }0 w (EtU = T (29) 
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and 



(£ T2 ) «(£ T2 ) = T 2 . (30) 

The above Lanczos-based methods for evaluating 4>\L<j)i are significantly more efficient than Chebychev-based 
methods used previously For a given level of accuracy in the present applications, they are between 3 and 5 

times more economical in the number of matrix multiplications required. Typically we achieve six figure convergence 
of the quadrature with 70 Lanczos steps on a matrix with a condition number (A max /A m i n ) of order 10 4 or 10 5 . 



Our goal was to achieve variance with respect to 4> (see (23)) which was one order of magnitude less than that with 



respect to the physical (gauge) distribution. We found that = 80 was a suitably conservative number of noise 
vectors to use. 

For estimating 5T, we note that 



E ST = E t , - E t = 4(L' ~ L)<fn ■ (31) 

r i=l 

Thus we can achieve variance 

a 2 (E ST ) = ^Tr((L'-L) 2 ) (32) 

if we use 

Est = Ej" — Et (33) 

provided we have employed the same set of noise vectors. This is simple to arrange. 
In fact, one could use stochastic estimators of the form (^oj) to estimate directly 

^ = ^Tr[l - Nr 1 ] = ^ReTr[l - M(M^M)- 1 ] . (34) 

An estimator for Tr[iW (M* M) -1 ] is obtained by setting 

v = 4>, u = M^<f> (35) 

where <j> is a suitable (e.g. Gaussian) noise vector. Then 

u^M^M^v = v^M{M^M)~\ (36) 

is an unbiased estimator of Tr[M(M^M) _1 ] as required. For this unsymmetric case [u ^ v), two Lanczos systems 
must be used and a subtraction performed Q. For the present analysis we have used the symmetric formalism as 
described above. Further analysis and discussion of these and related stochastic estimators is presented in fll|| . 
In the next section we report results of some tests of the matching procedures using the above algorithms. 



IV. NUMERICAL TESTS OF MATCHING 



A. Work estimates 



Having generated an ensemble of decorrelated configurations, one can consider estimating the numerical derivatives 
required for matching (see (|l5|)) in one of two ways: 

1. conduct two further simulations at neighbouring points in parameter space and make (uncorrelated) measure- 
ments of < F > on each of these 

2. apply the stochastic trace log techniques of the previous section to the existing ensemble. 
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One can estimate the relative amount of work involved in these two approaches. Suppose we seek to achieve an 
absolute error of e on a measurement of 5 < F > by each method and that the variances of < F > and of < 8TF > 
are a F and a 2 F respectively. The ratio of work required for the two approaches is then 

^2 = nr ajp W T , , 

Wi 4 ' a% ' W H mc 1 ' 

where Wt is the work done in a stochastic estimate of trace log on one configuration, tit is the number of trace logs 
required (either 2 or 3) and Wjjmc is the work done in generating one decorrelated configuration by Hybrid Monte 
Carlo (HMC). In turn, we estimate 

W T _ iV iV Lanc A 



Wumc N stc pN so i vc tac 

where the various parameters are associated with HMC simulation and the stochastic estimation of the determinant. 
These are defined in Table || which also shows typical values from the analysis presented below. With the numbers 
shown, the ratio Wt/Whmc is about 1/80. In Table |h| we show sample values of the variances up and a 2 F taken 
from the present study. Putting all this together we estimate 

^« 0.045. (39) 

Although the variance of < STF > can be quite large on large lattices, it is proportional to (Sk) 2 , or whatever the 
relevant difference parameter is. It is therefore possible, in principle, to obtain acceptable accuracy for the relevant 
derivative by method (2) with significantly less work. The above example suggests this is less than 5% of that required 
to obtain comparable accuracy by making additional simulations. One can easily refine the above treatment to take 
account of the work involved in equilibration. Of course, this makes method (2) seem even more attractive. 



B. Fixed plaquette curve 

First, we present results from matching the average plaquette < P >. This provides a check on some basic features 
of the procedures: the first order truncation of the cumulant expansion (||) and the Lanczos-based noisy estimator 
algorithm ( p6j ) . The initial reference point in the (3, k plane was taken as 

(0o, K ) = (5.2,0.1340) (40) 

where a sequence of well-equilibrated configurations had been generated by hybrid Monte Carlo on a 8 3 24 lattice 
with the non-perturbatively improved action for 2 flavours 0. At (3 = 5.2, the relevant improvement coefficient is 
c sw = 2.0171 | JTc| ] . It is estimated from preliminary spectroscopy measurements on larger lattices that K cr n is around 
0.1365 at this (3 and c sw - The values of k considered in this analysis are typical of those being used in production 
simulations. Plaquette measurements from all trajectories within this data sample showed autocorrelation times less 
than 20 trajectories. Longer runs showed an autocorrelation time for the plaquette of roughly double this. Bootstrap 
with binning was used to estimate the errors on all quan titie s. 



Using the algorithm and choices described in section III, we made stochastic estimates of T = TrmAf^M on a 



sequence of 40 configurations separated by 20 trajectories. In order to make use of ( |13| ) and (15) this must be done 
using a minimum of tit — 3 parameter sets: 

{c sw (A)),ko}, {c sw (/3 ),ko + 5k}, {c sw (I3 q + 5(3), k } . (41) 

Recall that for the unimproved Wilson action, there is no need to account for the additional (3 dependence in T which 
enters via c sw ((3). 

Having selected a change in bare quark mass Sk (for example Sk = —0.0005), we then use ( p"5| ) with F = P to 
estimate d/3/d.K. This yields a first estimate of the change 5(3 required to maintain a fixed value of < P >. One 
can then use (|l^) to verify this estimate of 5(3 by making further stochastic estimates of T at the final parameter set 
{c sw (Po + 5(3), k + 5k}. In all cases studied, this last verification step has been well satisfied and so one can in fact 
identify the matching curve directly from the two partial derivatives as proposed in the previous section. Results are 



shown in Table [II 



We have then generated further dynamical fcrmion configurations at the matched point 
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(/3o + SP, k + Sk) = (5.220, 0.1335) 



(42) 



and accumulated a similar ensemble of configurations for subsequent measurement. The corresponding value of c sw is 
1.9936. The plaquette measurements were made using relatively high statistics yielding the statistical errors shown. 
In order to make a proper comparison, one should fold into the error on < P > at (5.220, 1.9936, 0.1335) that due to 
the uncertainty in (3 (±0.0010). This would feed through to an additional uncertainty in < P > of ±0.0004. Thus 
the matching test is very well satisfied for the plaquette. 

The matching prediction done in reverse, back from n = .1335 to .1340, is also seen to be well satisfied. From Ta- 
ble |r| we see that (5.220,0.1335) is expected to match with (5.205(15) TU340) in good agreement with (5.2,0.1340) 



from where the matching estimates were originally made. Also in Table III, we show the estimated j3 shift correspond- 
ing to a further change of —0.0005 in n. Steps of this kind allow one to set up a grid of points from which which one 
can then deduce the fixed plaquette curve in the (3, k pla ne. 

For the sake of completeness, we also show in Table III the estimated shift required to match with quenched 
measurements of < P > (i. e. at infinite k). Independent gauge simulations at (3 = 5.61 show a good match of < P > 
when one takes into account the additional uncertainty in < P > of around ±0.009 which would feed through from the 
error of ±0.03 on estimating (3. Of course, for such large shifts the first order approximation may not be sufficiently 
accurate. We have calculated the second order approximation to 5(3 || for this quantity (±0.52(10)) but the statistical 
error is such that one cannot reliably discriminate it from the first order result (±0.41(3)) with the present level of 
statistics. 

A further plaquette matching test is given in Tabic |ll]. In this example, the reference point for the 'constant' 
plaquette curve was (5.2, 0.1330). Again, direct simulation showed that the matching was accurate and self-consistent. 
The matched points on this curve (< P >= 0.5197(3)), have been used to conduct tests of parallel tempering as 
described in section |v|. 



C. Full action matching 

We have noted in section || that some of the correlations required to identify curves of constant action ( |l8| ) are not 
directly calculable by the techniques of section III. Those of the form 

< T(cw, k)T{c' sw , k') > (43) 

require some care when setting up unbiased estimators. In particular, 

Ett' = Et^t' (44) 

is not unbiased. For an unbiased estimator, one requires something like that used for T 2 , i. e. E T 2. Unfortunately it 
is not so easy to evaluate the analogue of (^5|) via the above Lanczos methods. However, provided the variance of 
Et with respect to noise (^3|) is indeed small compared to that with respect to gauge fluctuations, the estimator ( (44| ) 
provides a useful approximation. Using this approximation, we have measured the shift 8 [3 corresponding to a shift Sk 



at fixed action to compare with that at fixed plaquette. For the first test shown in Table [II, we find 5(3 = 0.0200(10) 
(statistical error only) consistent with the value 0.0199(10) found for the fixed plaquette curve. In each case studied, 
we have found such consistency. 

We conclude that the fixed action and fixed plaquette curves are not significantly different at this order. 



D. Gauge invariant loops 

We are especially interested in matching those observables which are more sensitive to long range physics. We have 
repeated the above plaquette matching analysis using a variety of Wilson loops. The 16 loops used consist of 4 basic 
shapes realised in 4 different 'magnifications' (xl, x2, x3 and x4). The 4 basic shapes used were those specified in 
terms of link steps by the operations shown in Table [TV] and rotations thereof. 

For each loop L, we have evaluated the shift 5(3l required to hold the Wilson loop < Wl > constant under a change 
Sk (= —0.0005). The results of this are shown in Figure 0. The values of 5f3 L are similar to each other and to that for 
the average plaquette measured above. There is not much evidence of a shift increasing with the loop size, although 
the x2 loops show a higher trend than the plaquette value. One sees little evidence of mis-matching above the level 
of a standard deviation. 
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If this result (same 80) was reproduced for all gauge invariant loops and all linear combinations thereof, we would 
conclude that the static potential, ro and hence the lattice spacing, would be identical at the matched points ( J40| ) 
and ([42|). This would realise our initial objective of defining curves of constant effective volume. However, it cannot 
be that all ftxed-F curves emanating from a finite reference point (0q, kq) coincide. Moving from this point, in one 
direction we approach the quenched limit and in the other the chiral limit. We expect that different observables will 
be more or less sensitive to the effects of quenching. For example, the mass of a vector meson will probably change 
by less than 10% as the chiral limit is reached in the full theory as compared to its quenched value. On the other 
hand, the string tension should change from the lattice equivalent of 440 Mev to zero, eventually. We study the static 
potential in the next section. 



E. Potential and ro 

We have used the methods of Jl2| to measure the potential on each of the main ensembles studied. Since these 
are on 8 3 x 24 lattices, there are strong finite size effects present in V(r) and ro at the parameter values of interest. 
However, for the present purpose this is of little consequence. In the variational methods of jl2| one constructs 
'fuzzed' loops from a variety of spatial paths and employs transfer matrix methods to extract energy eigenvalues. The 
potential values were estimated by taking weighted averages of the effective masses at large time. We took care to 
use the same procedures on all ensembles. Errors were estimated by bootstrap. 

Figure || shows the static potential at the matched points p0| ) and ( |42| ) . The values are in good agreement at short 
distances but show a systematic divergence at larger separations. Figure || shows more clearly the difference between 
the potentials 1/(5.2,0.1340) — V(5.220, 0.1335). There appears to be a systematically increasing difference at larger 
distance. For comparison, the figure also shows V(5. 2, 0.1340) — 1^(5.2,0.1335) where there has been a shift in k but 
no compensating change in There is clear disagreement at all distances, as expected. 

The remaining set of points in Figure || shows the prediction for (5.220,0.1335) from the reference ensemble at 
(5.2,0.13400) using (0) to first order. Within the large statistical errors, the predicted difference is indeed compatible 
with zero. This demonstrates that where matching has been done only approximately (in this case with < P >) an 
observable can still be reliably estimated at another nearby point of interest. 

From the comparison of the directly simulated points, we conclude that matching the plaquette is not equivalent 
to matching the long range potential. Note that no corrections have been made for lattice artifacts at short distances 
but these are expected to be similar in each case. 

From the potential measurements represented in Figure § we can extract corresponding values of ro [Q and, hence, 
lattice spacing using ro(phys) = 0.49 fm. These are shown in Table |v|. In finding ro, we have calculated interpolated 
values of the static force via Newtonian 5-point interpolation of the potential. As usual, errors were computed via 
bootstrap. The same procedures were used for each ensemble. 

From Table we see that the lattice volume at the matched points is similar although that at the heavier quark 
mass is perhaps one standard deviation larger. This is a reflection of the observation made above that the potentials 
diverge at large separations due to their different slopes. We have attempted to estimate 8(3 required to match the 
mean values of the lattice spacings. In the last three rows of Table we show the lattice spacing predicted from 
the reference ensemble via ([|) for (3 shifts of +0.02, +0.03 and +0.04 to go with the k shift of —0.0005. The spacing 
predicted for 8(3 = 0.02 agrees well with that measured by direct simulation (previous line). From this value and 
those corresponding to 8(3 = 0.03 and 0.04 we estimate that the optimal matching for the lattice spacing (as opposed 
to the plaquette) would be 8(3 = 0.032(10). This value is compared with those corresponding to the various Wilson 
loops in Figure [l]. 

We conclude that the constant lattice volume curve (as defined by ro) may indeed differ from that corresponding 
to constant plaquette value. With only 40 configurations, the evidence is of marginal statistical significance. As a 
cross-check we measured r directly from a simulation at (5.232, .1335). See Table |y|. The result was compatible with 
the prediction and with the matched ensemble from which the prediction was made. Again, the statistical significance 
is not high. 

We have also measured ro on quenched configurations at (3 — 5.61, the point which is predicted and verified to have 
matched values of < P >, and at = 5.85, close to the point where ro is expected to match. Results are shown in 



Table [II. As expected, the lattice spacing does not match at 5.61 but is close to matching at 5.85. The shift required 



to match ro is some 60% larger than that required to match < P >. 
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F. Hadron correlators 



Finally in this section, we present results from matching lattice pion correlators. This test is made more practicable 
by recent advances in measuring hadron correlators with good statistical precision on a single gauge configuration fl3|| . 
We have made measurements of the local pion correlator C n (t) on the reference ensemble and calculated the corre- 
sponding 8(3 shifts. That is, for each value of t, we estimate the shift 8(3 required to keep constant for the test 



k shift. Results are presented Table VI where they are compared with examples of other observables. At short time 
separations (0 and 1) we find values compatible with that for the plaquette. At large values of t the results are 
overwhelmed by noise and we are unable to draw conclusions. However there are indications that for increasing t the 
shift required in (3 is also increasing. See, for example the correlator for t = 2 and corresponding effective mass values 
which are also shown in the table: 

mf(t) =-]n[C v (t + l)/C*(t)] . (45) 

Clearly one requires greater statistics to confirm these trends, but they are consistent with those inferred from the tq 
results presented above. A compilation of shifts 8/3 for sample loops, CV(i) and r are shown in Figure [|. 
Note that the effective mass values involved in this test are very far from those of a physical pion. 



V. FURTHER APPLICATIONS 



The above matching technology has a variety of potential uses including the following. 



A. Parallel tempering 



Parallel tempering (PT) is an improved Monte-Carlo method originally proposed by Hukushima et al 1 14 to improve 
simulations of spin glasses. It was further discussed by Marinari et al in |fl5f and fl^j who suggested its use in Lattice 
QCD. Recently Boyd |l7[] applied the technique to lattice QCD with staggered fermions and found evidence that 
Parallel Tempering did help decorrelate long distance observables. 

Parallel tempering essentially consists of running several independent simulations in parallel and with different 
parameters. Each such simulation produces a set of configurations which is distributed according to the probability 
distribution dictated by the simulation action and parameters. The PT algorithm exploits the fact that these distri- 
butions may have an overlap and occasionally attempts to swap configurations between ensembles. Acceptance of the 
swap is controlled by a Metropolis style acceptance step. 

This is the same situation described by the matching criterion M3 in section ||. From another viewpoint, the 
distance between the actions in matching criterion M2 can be related to the acceptance rate of the swaps in a parallel 
tempering algorithm. 

In frij| , tempering was carried out in the quark mass only. All the ensembles had the same value of (3. Furthermore, 
the quark masses had to be spaced quite closely together to obtain a reasonable swap acceptance rate. With the 
matching technology presented in this paper, it is possible to temper in both (3 and the quark mass. This may allow 
one to perform PT along a curve of approximately constant volume, and at a such a separation between ensembles 
that one might use some of the tempering ensembles to perform chiral extrapolations. Alternatively, one might be 
able to simulate with ensembles suitably chosen to improve the decorrelation properties of the system as a whole. A 
detailed investigation into PT using the matching technology is being conducted by us and the full results will be 
reported elsewhere fL8[. 



B. Approximate algorithms 

In H we demonstrated how the parameters of approximate algorithms could be tuned according to the criteria 
Ml, M2 or M3 described above. In subsequent tests @, we showed that approximate algorithms based on a few 
Wilson loops only, were unlikely to produce a very accurate approximation, at least in the sense of M3 where the 
approximate action acts as the guide within an exact algorithm. The variance of the difference between the two 
extensive quantities remains unacceptably large on lattices of a useful size. It is, however, still of considerable interest 
to design approximate or model actions which improve on the quenched approximation by encapsulating at least some 
of the additional physics implied by dynamical quarks. 
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An alternative route to an exact algorithm might be to make use of the Lanczos quadrature approximation of 
section ^ for part of the effective action and gauge invariant loops as above for the remainder. The trial configurations 
would be generated by the loop part of the action and an accept/reject step based on the Lanczos part. The technology 
of section [□] can be used to tune the loop part together with the Lanczos part to match the exact action. As a simple 
example of this approach, consider an approximate action defined such that (c.f. (0)) 

S app = -P'Wa - T(A Lanc ) (46) 



where T(Ni janc ) is the approximation to Trln(AfTM) as described in section [II but using only ALanc Lanczos itera- 
tions. The loop part of the effective action is just a single plaquette in this example. In a standard Metropolis update 
this action would only be viable if A^anc was considerably smaller than the typical values (90) required to estimate 
the true action, and sufficient account was taken of the short range fluctuations by having a properly tuned gauge 
loop part. This approximation is similar in spirit to that advocated in [Is|j20[| where it is argued that a truncated sum 
of low- lying eigenvalues can reproduce the gross behaviour of the fermion determinant. In the present scheme, we are 
able to obtain a particularly efficient approximation to the trace log by using the optimal weighting determined by 
the Gaussian quadrature rule. We have conducted preliminary tests of these ideas by measuring the shift 6(3 — (3' — (3 
required to compensate for a truncation to A^anc Lanczos iterations. In a simple test which matched the average 
plaquette, the shift in (3 was reduced by a factor of 8 in changing from A^anc = (quenched) to A^anc = 4, for 
example. The correspond residual variance (after matching) also dropped by a factor of around 8. For increasing 
A^Lanc, the shift rapidly becomes compatible with zero at the level of statistical accuracy implied by the number of 
noise vectors used (80). This demonstrates that the Lanczos quadrature approach gives a very efficient estimator for 
the trace log. It will be worth exploring whether the long range modes described by such an approximation can be 
combined with a suitably tuned gauge loop action describing short range modes so as to achieve a practical exact 
algorithm. 



VI. CONCLUSIONS 



We have proposed a strategy for dynamical quark simulations in which the effective lattice volume is held fixed 
while the effects of progressively lighter sea quarks are investigated. Possible algorithms for accomplishing this have 
been presented and the results of tests discussed. In particular, we have presented results using an efficient stochastic 
estimator of the fermion determinant and quantities related to it. These include estimates of the constant lattice 
spacing curves at relevant points in the (3, k plane for lattice QCD using a non-perturbatively improved action. 
We have demonstrated that the work involved in determining such curves via our differential stochastic methods is 
considerably less than that required to establish them by direct simulation. Further applications of these tcchniquess 
have been discussed. 
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FIG. 1. Predictions for the shift 8/3 from reference point (5.2,0.1340) obtained using the 16 sample loops described in the 
text. The diamond point is the corresponding value deduced from ro- The closed (open) points correspond to a shift of 
8k = -0.0005(+0.0005) from ensembles at (5.2,0.1340) and (5.220,0.1335) respectively. 
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FIG. 3. Differences in the static potential measured from the reference point (5.2, 0.1340). Solid and open circles correspond 
to (5.220,0.1335) and (5.2,0.1335) respectiveley. Open diamonds corespond to the prediction for (5.220,0.1335) using (^). 
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FIG. 4. Predictions for the shift 8(3 from reference point (5.2,0.1340) required to keep the specified observables constant 
when k is changed by 8k = —0.0005. Circles correspond to Wilson loops (1 x 1, 1 x 2, 2 x 4); diamonds to C,r(i) (t = 0, 1, 2); 
square to r . 
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TABLE I. Parameters association with HMC simulation and with the stochastic estimation of the determinant. 



Parameter 


Description 


Value 


iVstep 


number of HMC steps per trajectory 


50 




number of sweeps in the HMC solver (e.g.BiCGStab) 


300 


A 


HMC trajectory acceptance 


0.75 


1~AC 


autocorrelation time in trajectories 


30 



Nj, number of noise vectors (section |II]|) 80 

iVLanc number of Lanczos iterations 90 



TABLE II. Sample variances usued to estimate the relative work involved in measuring differences directly and by the 
techniques proposed in this paper. These correspond to a shift Sk = —.0005 and the reference data described in the text. 

TTi „2 2 2 /2 

£ £> ^SF vsf/vf 

P (average plaquette) 5.2 x 10~ b 2.9 x 10~ b 5J5 

r 0.78 2.7 3.5 



TABLE III. Matching the average plaquette. 



(A), Csw, Ko) 


< P > 


traj 


K 


5 k 


8(3 


config 


(5.200,2.0171,0.1340) 


0.5286(3) 


6000 


0.1335 


-0.0005 


+0.0199(10) 


40 













+0.41(3) 


40 


(5.220,1.9936,0.1335) 


0.5290(3) 


6000 


0.1340 


+0.0005 


-0.0195(15) 


40 








0.1330 


-0.0005 


+0.0162(23) 


40 













+0.40(4) 


40 


(5.61,0) 


0.5275(3) 


1000 sweeps 










(5.200,2.0171,0.1330) 
(5.215, 1.9994, 0.1325) 


0.5197(3) 
0.5207(6) 


6000 
2000 


0.1325 



+0.015(3) 
+0.368(16) 




40 
40 



TABLE IV. Construction of sample gauge invariant loops from specified link steps. A step +2 means a link along lattice 
direction 2 while —1 means a link along lattice axis 1 in the negative direction. 



Loop No. of links Link steps 

~ 4 (+l,+2,-l,-2) 

2 6 (+l,+l,+2,-l,-l,-2) 

3 6 (+l,+2,+3,-2,-l,-3) 

4 6 (+l,+2,+3,-l,-2,-3) 



TABLE V. Values of ro and lattice spacing deduced from the static potential. The last two rows are values predicted by 
from the reference ensemble (5.2,0.13400). 



03, k) 


ro/a 


a fm 




(5.200,0.1340) 


3.87(17) 


0.127(6) 


direct simulation 


(5.220,0.1335) 


3.69(11) 


0.133(4) 


by HMC 


(5.232,0.1335) 


3.76(13) 


0.130(4) 




(5.61,0) 


2.33(2) 


0.211(2) 




(5.220,0.1335) 


3.67(18) 


0.134(6) 


pred. 


(5.230,0.1335) 


3.80(22) 


0.129(9) 


from (5.2^0.13400) 


(5.240,0.1335) 


4.06(33) 


0.121(12) 


using (@) 
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TABLE VI. Summary of the shifts <5/3 required to maintain the specified quantities constant under a change 8k = —0.0005 
from the reference point (f3 , /to) = (5.200, 0.13400). 



Measurement 


Value at [po,Ko) 


OP 


< r > 




n m no 1 1 n\ 

u.uiyy^iuj 


< Wlx2 > 


0.3161(4) 


0.0207(12) 


< W 2X 4 > 


0.0312(2) 


0.0213(15) 


ro /a 


3.87(17) 


0.0302(10) 


CV(0) 


2.132(2) 


0.0230(18) 


CV(1) 


0.2782(8) 


0.026(3) 


CV(2) 


0.0792(5) 


0.030(5) 


C„(3) 


0.0298(3) 


0.035(80) 


< ff (l) 


2.037(2) 


0.030(5) 


<ff(2) 


1.256(4) 


0.035(28) 
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